7-20-22 stable-gfp pup1-rfp decay

Only 20min timelapse because we concluded that we wont be able to decouple photobleaching from proteasomal decay

this RMD is to generate the df for decay rate estimation for the DMSO control

library(data.table)
library(tidyverse)
library(reshape2)
library(gganimate)
library(ggridges)
library(ggpubr)

source("~/R-scripts/R_functions/master-2-3-22.R")

#getting the information for imaris files ##files to read from Imaris

FilesToRead <- c("C0_Area", 
            "C0_Ellipsoid_Axis_Length_B",
            "C0_Ellipsoid_Axis_Length_C",
            "C0_Intensity_Center_Ch=3",
            "C0_Intensity_Mean_Ch=3", 
            "C0_Intensity_Median_Ch=3",
            "C0_Intensity_Sum_Ch=3",
            "C0_Number_of_Voxels", 
            "C0_Number_of_Triangles",
            "C0_Sphericity", 
            "C0_Volume.csv", 
            "C0_Position.csv" )

##actual reading and making them into df

####MG131
#Getting the list of all the samples in the experiment 
gfp_imaris <- dir("~/Imaris-data/2022/pup1_rfp/7-20-22-pup1/stable/imaris/gfp", full.names = TRUE)

#getting the information for imageJ image info files
#getting all the image info txt for all the images in the experiment. 
list.image.info <- list.files("~/Imaris-data/2022/pup1_rfp/7-20-22-pup1/stable/image_info_imageJ" , full.names = TRUE)

#Extracting information from the imageJ csv into a df with actual time differences between each images 
real.time.interval.df <- get.image.info.OnlyGFP(list.image.info = list.image.info)      

#Extracting information from the imaris csv into a df with information on area, time, gfp intensity, mcherry intensity, (mean, median and sum), volume, trackID and IDs
list.of.files.gfp <- get.imaris.info.onlyGFP( dir.names = gfp_imaris, 
                                              files_to_read = FilesToRead ,
                                              real.time.interval.df = real.time.interval.df)

list.of.files.gfp$`20min` %>% 
  split(.$field)

#Assigning the sample by to each field. There are two fields of autofluorescecne in this experiment. The neg.singal stands for the autofluorescence sample.

list.of.files.gfp[["20min"]] <- list.of.files.gfp[["20min"]] %>% 
  mutate(sample = ifelse(exp.field %in% c("20min_s1", "20min_s2"), "neg.signal", "pos.signal"))

#getting the imaris info of DAPI and pup1-rfp ##files to read (these will have rfp and dapi) ##Use the cells_surface data from Imaris #files to read

#reading the cell surface imaris files

# dapi_pup1_imaris_cells <- dir("~/Imaris-data/2022/proteasome_inhibition/2-23-22-PI/1uM/imaris/surface/cells", full.names = TRUE) 
# 
# list.of.files.dapi.pup1 <- get.imaris.info.dapi_pup1(dir.names = dapi_pup1_imaris_cells, 
#                                                      files_to_read = FilesToReadDAPIRFP)
# 
# unique(list.of.files.dapi.pup1$timepoint)

#getting the pup1-rfp puncta info ##use the pup1_surface info from Imaris

FilesToReadDAPIRFP <- c("C0_Area", 
            "C0_BoundingBoxOO_Length_B",
            "C0_BoundingBoxOO_Length_C",
            "C0_Ellipsoid_Axis_Length_B",
            "C0_Ellipsoid_Axis_Length_C",
            "C0_Intensity_Center_Ch=4",
            "C0_Intensity_Mean_Ch=4", 
            "C0_Intensity_Median_Ch=4",
            "C0_Intensity_Sum_Ch=4",
            "C0_Intensity_Center_Ch=5",
            "C0_Intensity_Mean_Ch=5", 
            "C0_Intensity_Median_Ch=5",
            "C0_Intensity_Sum_Ch=5",
            "C0_Number_of_Voxels", 
            "C0_Number_of_Triangles",
            "C0_Sphericity", 
            "C0_Volume.csv", 
            "C0_Position_X.csv",
            "C0_Position_Y.csv",
            "C0_Intensity_Max_Ch=6",
            "C0_Intensity_Center_Ch=6")

#first get the info of the cell surface then match each cell surface with the pup1 puncta. do not include the MRG cells as they have no pup1-rfp puncta mask

dapi_pup1_imaris_cells <- dir("~/Imaris-data/2022/pup1_rfp/7-20-22-pup1/stable/imaris/surface/cells", full.names = TRUE) 

list.of.files.dapi.pup1 <- get.imaris.info.dapi_pup1.surface(dir.names = dapi_pup1_imaris_cells, files_to_read = FilesToReadDAPIRFP) #updated this function to not include trackIDs

list.of.files.dapi.pup1 <- list.of.files.dapi.pup1 %>%
    mutate(timepoint = ifelse(timepoint == 2,31,timepoint))
list.of.files.dapi.pup1 %>% split(.$timepoint)

#Getting the pup1 puncta surface info

pup1_surface_imaris <- dir("~/Imaris-data/2022/pup1_rfp/7-20-22-pup1/stable/imaris/surface/pup1", full.names = TRUE) 

pup1_puncta <- get.imaris.info.dapi_pup1.surface(dir.names = pup1_surface_imaris[3:8], 
                                                 files_to_read = FilesToReadDAPIRFP) 
pup1_puncta <- pup1_puncta %>% 
  mutate(timepoint = ifelse(timepoint == 2,31,timepoint))

pup1_puncta %>% split(.$timepoint)

#combine the cell surface info with the pup1 puncta info.

#this should have 11159 no. of cells 
pup1.dapi.cell.puncta <- list.of.files.dapi.pup1 %>% #df with entire cell segmented 
  select(
    area,
    pos.x,
    pos.y,
    timepoint,
    all.mask.int.max,
    all.mask.int.center,
    rfp.int.mean,
    rfp.int.median,
    rfp.int.sum,
    dapi.int.mean,
    dapi.int.median,
    dapi.int.sum,
    experiment,
    field
  ) %>% 
  left_join(
    .,
    pup1_puncta %>%        #df with the pup1 puncta segmented 
      select(
        rfp.int.center,
        rfp.int.mean,
        rfp.int.median,
        rfp.int.sum,
        dapi.int.mean,
        dapi.int.sum,
        dapi.int.median,
        all.mask.int.max,
        all.mask.int.center,
        area,
        timepoint,
        experiment,
        field,
        no.of.voxels
      ) %>% rename("area.puncta" = "area",
                   "rfp.center.puncta" = "rfp.int.center",
                   "rfp.mean.puncta" = "rfp.int.mean",
                   "rfp.sum.puncta" = "rfp.int.sum",
                   "rfp.median.puncta" = "rfp.int.median",
                   "dapi.mean.puncta" = "dapi.int.mean",
                   "dapi.sum.puncta" = "dapi.int.sum",
                   "dapi.median.puncta" = "dapi.int.median",
                   "no.of.voxels.puncta" = "no.of.voxels"),
     by = c("all.mask.int.max" , 
            "experiment" , 
            "field" , 
            "timepoint" , 
            "all.mask.int.center")
  )


  
#this has 6514 cells. some cells have multiple puncta. maybe just drop the cells. 

#combine the cell surface + pup1 puncta info of the pup1 expressing cells with all the cell surface information. There will be NAs in the columns with puncta in their names for the neg cells

# pup1.dapi.all <- list.of.files.dapi.pup1 %>% 
#   left_join(.,pup1.dapi.cell.puncta %>% 
#               select(.,area, 
#                      pos.x, 
#                      pos.y,
#                      timepoint,
#                      experiment, 
#                      field, 
#                      rfp.center.puncta, rfp.mean.puncta, rfp.median.puncta, rfp.sum.puncta,
#                      dapi.mean.puncta,dapi.sum.puncta, dapi.median.puncta,
#                      area.puncta, no.of.voxels.puncta), 
#             by = c("area","pos.x","pos.y","timepoint", "experiment", "field"))
# 
# pup1.dapi.all

#dont need to do this. 

Combining the GFP intensities with the rfp puncta info at t = 1

cell_atributes <- list.of.files.gfp$`20min` %>% 
  filter(timepoint %in% c(1,31)) %>% 
  left_join(.,pup1.dapi.cell.puncta %>% select(-all.mask.int.max,
                                               -all.mask.int.center) ,
            by = c("area",
                   "timepoint",
                   "pos.x",
                   "pos.y",
                   "experiment",
                   "field"))

cell_atributes %>% 
  # filter(unique.trackID == "1000001774_20min_s3")
    dplyr::group_by(field, unique.trackID, timepoint) %>%
    dplyr::summarise(n = dplyr::n(), .groups = "drop") %>%
    dplyr::filter(n > 1L)

#6440 cells. There are 6 cells with multiple pup1 puncta. maybe drop these cells just for simplicity purposes. 

#Background subtraction

#background info files for gfp and mcherry and DAPI: the background folder has two directories: gfp bg and mcherry+dapi bg
bg <- dir(path = "~/Imaris-data/2022/pup1_rfp/7-20-22-pup1/stable/background" , full.names = TRUE)

#Get the information for the background intensity for every image for 
#rfp and DAPI
rfp.dapi.bg.means <- bg.intensity.rfp.dapi(bg.files = list.files(bg[1] , full.names = TRUE))
rfp.dapi.bg.means

$`20min`

#for GFP
gfp.bg.mean <- bg.intensity.gfp(bg.files = list.files(bg[2] , full.names = TRUE))

gfp.bg.mean$`20min` %>% arrange(field,timepoint)

#gfp background for s4, s5 and s6 samples:

#actually subtract the background intensity from the gfp and mcherry intensity of the cells

#first fun. sub.bg.intensity subtracts the bg intensity from each cell. 

gfp.bg.sub <- sub.bg.intensity.gfp(list.of.files = list.of.files.gfp, 
                                   bg.means = gfp.bg.mean$`20min`  )
gfp.bg.sub

$`20min`

rfp.dapi.bg.sub <- sub.bg.intensity.dapi.rfp(list.of.files = cell_atributes, 
                                             bg.means = rfp.dapi.bg.means$`20min` )

gfp.bg.sub$`20min` %>% 
  filter(gfp.mean.bg.sub > 0,
         ) %>% 
  group_by(field, timepoint) %>% 
  tally() 

rfp.dapi.bg.sub$`20min` %>% split(.$field)

$s1

$s2

$s3

$s4

$s5

$s6

$s7

$s8

NA

#combine the background subtracted gfp dataframe with the bg subtracted rfp/dapi total and puncta dfs

gfp.rfp.dapi.bg.sub <- gfp.bg.sub$`20min` %>% 
  filter(timepoint %in% c(1,31)) %>% 
  select(unique.trackID,
         gfp.mean.bg.sub,
         gfp.sum.bg.sub,
         avg.gfp.bg,
         area,
         pos.x,
         pos.y) %>% 
  left_join(.,rfp.dapi.bg.sub$`20min`, by = c("unique.trackID", "area","pos.x","pos.y")) %>% split(.$experiment)

#filtering cells which have negative background subtracted intensities

#Uniqe trackIDs of cells above the background int for gfp, dapi and rfp 
#this function results in a list of two dfs. the 1st df = uniquetrackIDs with the experiment they're from.
listOfUniqtrackID_gfpfiltereddf <- trackID.abv.zero.gfp.rfp.dapi(
  bg_subtracted_dftp1 = gfp.rfp.dapi.bg.sub,
  list_gfp_all_tps = gfp.bg.sub)

#the df with gfp intensities for all timepoints which is filtered for cells gfp.mean.bg.sub > 0 
list.of.files.bg.sub.above.0.gfp <- listOfUniqtrackID_gfpfiltereddf[[2]] %>% 
  split(.$experiment)

#removing cells with negative background subtracted intensities from the dataframe with the cell info from timepoint 1 and the last timepoint

gfp_rfp_dapi_bgsubAbvZero_tp1 <- listOfUniqtrackID_gfpfiltereddf[[1]] %>% 
  left_join(.,bind_rows(gfp.rfp.dapi.bg.sub) , by = c("unique.trackID", "experiment"))

gfp_rfp_dapi_bgsubAbvZero_tp1 %>% 
  group_by(timepoint, exp.field) %>% 
  tally()

list.of.files.bg.sub.above.0.gfp$`20min`%>% 
  group_by(timepoint, exp.field) %>% tally()

#s3, s4, s5, s6, s7 and s8 need to fix the background subtraction

#density plot of the DAPI staining. Higher DAPI staining == dead cells

#at the first timepoint
gfp_rfp_dapi_bgsubAbvZero_tp1  %>% 
  ggplot(.,aes(x = dapi.mean.bg.sub, color = as.factor(timepoint)))+
  geom_density(aes(y=..scaled..))+
  facet_wrap(~field, scales = "free_x")+
  scale_x_log10()+
  theme_pubr(x.text.angle = 45)

#correlation of the dapi staining at timepint 1 and the last timepoint
# gfp_rfp_dapi_bgsubAbvZero_tp1  %>%
#   select(dapi.mean.bg.sub, timepoint, field, unique.trackID) %>%
#   filter(!(unique.trackID == "1000000259_20min_s8")) %>% 
#   pivot_wider(values_from = dapi.mean.bg.sub, names_from = timepoint) %>% 
#   ggplot(.,aes(x = `1`, y = `31`))+
#   geom_point()+
#   facet_wrap(~field, scales = "free_x")+
#   theme_pubr(x.text.angle = 45)

#the above code did not work because there are duplicate values for a given cell ID. 

#to find the duplicates
gfp_rfp_dapi_bgsubAbvZero_tp1  %>% 
  select(dapi.mean.bg.sub, timepoint, field, unique.trackID) %>%
    dplyr::group_by(field, unique.trackID, timepoint) %>%
    dplyr::summarise(n = dplyr::n(), .groups = "drop") %>%
    dplyr::filter(n > 1L)

#remove the duplicates and look at how dapi staining is changing from timepoint 1 to tp 31:
#correlation of the dapi staining at timepint 1 and the last timepoint
gfp_rfp_dapi_bgsubAbvZero_tp1  %>%
  select(dapi.mean.bg.sub, timepoint, field, unique.trackID) %>%
  filter(!(unique.trackID == "1000002071_20min_s8")) %>% 
  pivot_wider(values_from = dapi.mean.bg.sub, names_from = timepoint) %>% 
  ggplot(.,aes(x = `1`, y = `31`))+
  geom_point()+
  facet_wrap(~field, scales = "free_x")+
  theme_pubr(x.text.angle = 45)

Warning: Values from `dapi.mean.bg.sub` are not uniquely identified; output will contain list-cols.
* Use `values_fn = list` to suppress this warning.
* Use `values_fn = {summary_fun}` to summarise duplicates.
* Use the following dplyr code to identify duplicates.
  {data} %>%
    dplyr::group_by(field, unique.trackID, timepoint) %>%
    dplyr::summarise(n = dplyr::n(), .groups = "drop") %>%
    dplyr::filter(n > 1L)
Error: Discrete value supplied to continuous scale

#get the threshold of dapi above which the cells will be deemed dead!!

#removing cells with dapi intensity greater than the dapi threshold set

#getting the IDs of the cells remaining after removing the cells with higher than threshold dapi intensity. use this to filter cells from the df with mcherry and gfp timelapse information
#will remove dead cells from the first timepoint 
dapi.filtered.cellsTrackID <- bind_rows(deadcell_temp_df) %>% 
  filter(timepoint %in% c(1,31)) %>% 
  left_join(.,deadCell_boundry, by = "exp.field") %>% 
  filter(dapi.mean.bg.sub < dapi.threshold) %>% pull(unique.trackID)

#df with fluorescent intensities of gfp, rfp and dapi for live cells at tp =1 ant tp = 31.
liveCellsRFP.DAPI.GFP <- gfp_rfp_dapi_bgsubAbvZero_tp1 %>% 
  left_join(.,deadCell_boundry, by = "exp.field") %>% 
  filter(dapi.mean.bg.sub < dapi.threshold)


#df with dead cells removed this df has gfp intensities for all the timepoint which are filtered for dead cells 
new.gfpLiveCells <-  bind_rows(list.of.files.bg.sub.above.0.gfp) %>% 
  filter(unique.trackID %in% liveCellsRFP.DAPI.GFP$unique.trackID)

#saving the dataframe with the background subtracted, dead cells filtered intensities of gfps from the autofluor cells (MRG) and the gfp positive cells.

looking at the distribution of rfp and gfp in the live cells

#plotting line plots raw without any AF filtering

bind_rows(new.gfpLiveCells)%>%
  # filter(field %in% c("s7","s8")) %>% 
  ggplot(.,aes(timepoint, gfp.mean.bg.sub, group = unique.trackID ,color = sample))+
  geom_line(alpha = 0.2)+
  facet_wrap(~field, scales = "free_y")

# ggsave(filename="ln-mrg_vs_gfp.png", plot= last_plot() , path = "~/plots/7-2-20/plots_thesis_com_meet-2020/" , width = 13 , height = 8, bg = "transparent")

new.gfpLiveCells %>% 
  ggplot(.,aes(real.time.gfp, gfp.mean.bg.sub, group = unique.trackID ,color = sample))+
  geom_line(alpha = 0.2)+
  facet_wrap(~field, scales = "free_y")

  # annotation_logticks()

Fix background for s1 and s8: fixed background intensities of s1, s3, 4 till s8.

#Filtering cells based on autofluorescence

##getting the autofluorescence threshold for gfp and mcherry

#filtering based on pup1-rfp expression ##filtering cells which are above the gfp and rfp threshold In the filter.cells.below.autofluor.31tp.gfpV2() function, the function looks for the column named “threshold” and I need to decide which column will be threshold based on the overlap of the neg and positive cells. This is because there are cells which become dead/become in focus as you keep taking images which might suddenly increase the 95th quantile intensity of the negative cells.

#get a filtered df of cells passing the gfp threshold set above 
#this function below first removes cells with pup1-rfp less than the threshold set. 
#Then it removes the cells with less than AF threshold set at tp =1 
#After that it subtracts the af value from every cell for every timepoint, and if the subtracted value is less than 1 it makes it NA. 

cells.filtered.df.31tp <- filter.cells.below.autofluor.31tp.gfpV2(
  df_rfp_dapi_gfp = liveCellsRFP.DAPI.GFP,
  
  rfp_min_int = rfp.min.int %>% rename("quant_95" = "threshold_95", 
                                       "threshold_95" = "pup1_threshold"),
  
  gfp_live_cells_list = new.gfpLiveCells %>% split(.$experiment),
  
  gfp_min_int = gfp.min.int %>% 
    mutate(threshold = threshold_95 )
)

filtered.df.new.wo.af.31tp <- cells.filtered.df.31tp$filtered.df.new %>% 
  bind_rows() %>% 
  filter(sample == "pos.signal") %>% 
  split(.$experiment)

unique(filtered.df.new.wo.af.31tp$`20min`$timepoint)

 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31

remove the first two blurry images from s7

# filtered.df.new.wo.af.31tp$`20min` <- filtered.df.new.wo.af.31tp$`20min` %>% filter(field == "s7")
#   filter(!(field == "s7" & timepoint %in% c(1,2))) 

#making df for bulk regression

###This function creats the df for fitting the 24 models 
# list.df.model <- df.model(filtered.df.new = filtered.df.new.wo.af)
df.model.31tp.gfp <- df.model.gfp(filtered.df.new = filtered.df.new.wo.af.31tp)
df.model.31tp.gfp

####This function cleans up the df where you subtract the t=0 intensity with all the other intensities for gfp
final.data.31tp.gfp <- clean.df.modeling.gfp(data.gfp = df.model.31tp.gfp)

On 1/28/21 talked to premal about filtering: 1. Remove cells which have positive delta Intensity more than 2 times

#getting IDs of cells which have positive intensity values for less than 2 timepoints 
ids.int.2 <- final.data.31tp.gfp %>% 
  group_by(unique.trackID) %>% 
  tally(ln.gfp.dif > 0) %>% 
  filter(n<2) %>%            #keep cells which have less than two positive intensity values
  pull(unique.trackID)

#filtering cells which have trackID in the above vector  
#gfp
temp.posIntBlw2.gfp <- final.data.31tp.gfp %>% 
  filter(unique.trackID %in% ids.int.2) 

df_for_MechanisticModel <- temp.posIntBlw2.gfp  %>%
  left_join(
    .,
    bind_rows(filtered.df.new.wo.af.31tp) ,
    by = c(
      "unique.trackID",
      "timepoint" ,
      "experiment" ,
      "gfp.mean.bg.af.sub.new" ,
      "exp.field"
    )
  )

#save this dataframe for half life estimation

#saving the data

#SAVING THE DATAFRAME WHICH HAS ALL THE BACKGROUND INT, AUTOFLUOR INTENSITY (AT THE INTERSECTION, THE 80TH AND THE 95TH QUANTILE OF MRG AUTO GFP INTENSITY)